
#include "gsl_ext_function.h"
#include "gsl_ext_symbolic_function.h"

#include <gsl/gsl_errno.h>
#include <matheval.h>

gsl_ext_multimin_symbolic_function_f_evaluator *gsl_ext_multimin_symbolic_function_f_evaluator_alloc(char *expr)
{
  gsl_ext_multimin_symbolic_function_f_evaluator *e = malloc(sizeof(gsl_ext_multimin_symbolic_function_f_evaluator));

  e->f_evaluator = evaluator_create(expr);
  if(!e->f_evaluator)
  {
    free(e);
    GSL_ERROR_VAL("invalid symbolic function expression", GSL_EBADFUNC, 0);
  }

  evaluator_get_variables(e->f_evaluator, &e->var_names, &e->n);

  return e;
}

void gsl_ext_multimin_symbolic_function_f_evaluator_free(gsl_ext_multimin_symbolic_function_f_evaluator *e)
{
  evaluator_destroy(e->f_evaluator);
  free(e);
}

size_t 
gsl_ext_multimin_symbolic_function_f_evaluator_n(gsl_ext_multimin_symbolic_function_f_evaluator *e)
{
  return e->n;
}

gsl_ext_multimin_symbolic_function_df_evaluator *gsl_ext_multimin_symbolic_function_df_evaluator_alloc(gsl_ext_multimin_symbolic_function_f_evaluator *f_eval)
{
  int i;
	
  gsl_ext_multimin_symbolic_function_df_evaluator *e = malloc(sizeof(gsl_ext_multimin_symbolic_function_df_evaluator));
	
  e->f_evaluator = f_eval;
  e->df_evaluators = malloc(f_eval->n * sizeof(void *));
	
  for(i = 0; i < f_eval->n; i++)
    e->df_evaluators[i] = evaluator_derivative(f_eval->f_evaluator, f_eval->var_names[i]);
	
  return e;
}

void gsl_ext_multimin_symbolic_function_df_evaluator_free(gsl_ext_multimin_symbolic_function_df_evaluator *e)
{
  int i;
  
  for(i = 0; i < e->f_evaluator->n; i++)
    evaluator_destroy(e->df_evaluators[i]);
  
  free(e->df_evaluators);
  free(e);
}

size_t 
gsl_ext_multimin_symbolic_function_df_evaluator_n(gsl_ext_multimin_symbolic_function_df_evaluator *e)
{
  return e->f_evaluator->n;
}

gsl_ext_multimin_symbolic_function_d2f_evaluator *gsl_ext_multimin_symbolic_function_d2f_evaluator_alloc(gsl_ext_multimin_symbolic_function_df_evaluator *df_eval)
{
  int i, j;
  const int n = df_eval->f_evaluator->n;
  
  gsl_ext_multimin_symbolic_function_d2f_evaluator *e = malloc(sizeof(gsl_ext_multimin_symbolic_function_d2f_evaluator));
  
  e->df_evaluator = df_eval;
  e->d2f_evaluators = malloc(n * sizeof(void **));
  
  for(i = 0; i < n; i++)
    {
      e->d2f_evaluators[i] = malloc(n * sizeof(void *));
      for(j = 0; j <= i; j++)
	e->d2f_evaluators[i][j] = evaluator_derivative(df_eval->df_evaluators[i], df_eval->f_evaluator->var_names[j]);
    }
  
  return e;
} 

void gsl_ext_multimin_symbolic_function_d2f_evaluator_free(gsl_ext_multimin_symbolic_function_d2f_evaluator *e)
{
  int i, j;
  const int N = e->df_evaluator->f_evaluator->n;;
  
  for(i = 0; i < N; i++)
    for(j = 0; j <= i; j++)
      if(j <= i)
        evaluator_destroy(e->d2f_evaluators[i][j]);
  
  for(i = 0; i < N; i++)
    free(e->d2f_evaluators[i]);
  free(e->d2f_evaluators);
  free(e);
}

size_t 
gsl_ext_multimin_symbolic_function_d2f_evaluator_n(gsl_ext_multimin_symbolic_function_d2f_evaluator *e)
{
  return e->df_evaluator->f_evaluator->n;
}

double gsl_ext_multimin_symbolic_function_f(const gsl_vector *x, void *params)
{
  gsl_ext_multimin_function_params *p = (gsl_ext_multimin_function_params *)params;
  gsl_ext_multimin_symbolic_function_f_evaluator *eval = (gsl_ext_multimin_symbolic_function_f_evaluator *)p->f_evaluator;
	
  p->f_eval_counter++;
  
  return evaluator_evaluate(eval->f_evaluator, eval->n, eval->var_names, x->data);
}

void gsl_ext_multimin_symbolic_function_df(const gsl_vector *x, void *params, gsl_vector *g)
{
  int i;
  double gi;
  gsl_ext_multimin_function_params *p = (gsl_ext_multimin_function_params *)params;
  gsl_ext_multimin_symbolic_function_df_evaluator *eval = (gsl_ext_multimin_symbolic_function_df_evaluator *)p->df_evaluator;
	
  for(i = 0; i < eval->f_evaluator->n; i++)
    {
      gi = evaluator_evaluate(eval->df_evaluators[i], eval->f_evaluator->n, eval->f_evaluator->var_names, x->data);
      gsl_vector_set(g, i, gi);
    }
  
  p->df_eval_counter++;
}

void gsl_ext_multimin_symbolic_function_d2f(const gsl_vector *x, void *params, gsl_matrix *d2f)
{
  if(d2f->size1 != d2f->size2)
    GSL_ERROR("The Hessian matrix must be a square matrix.", GSL_EINVAL);
  if(x->size != d2f->size1)
    GSL_ERROR("The vector x is incompatible with the Hessian matrix.", GSL_EINVAL);
  
  const int N = d2f->size1;
  
  int i, j;
  double Hij;
  
  gsl_ext_multimin_function_params *p = 
    (gsl_ext_multimin_function_params *)params;
  gsl_ext_multimin_symbolic_function_d2f_evaluator *eval = 
    (gsl_ext_multimin_symbolic_function_d2f_evaluator *)p->d2f_evaluator;
  
  if(x->size != gsl_ext_multimin_symbolic_function_d2f_evaluator_n(eval))
    GSL_ERROR("The vector x is incompatible with the function.", GSL_EINVAL);
  
  /* NOTE: only the lower-triangular part is computed */
  for(i = 0; i < N; i++)
  {
    for(j = 0; j < N; j++)
    {
      if(j > i)
        continue;
      Hij = evaluator_evaluate(eval->d2f_evaluators[i][j], N, 
                               eval->df_evaluator->f_evaluator->var_names, 
                               x->data);
      gsl_matrix_set(d2f, i, j, Hij);
    }      
  }
  
  p->d2f_eval_counter++;
}
